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Abstract 

An optimized molecular model for ammonia, which is based on a previ- 
ous work of Kristof et al., Mol. Phys. 97 (1999) 1129-1137, is presented. 
Improvements are achieved by including data on geometry and electro- 
statics from ab initio quantum mechanical calculations in a first model. 
Afterwards the parameters of the Lennard- Jones potential, modeling dis- 
persive and repulsive interactions, are optimized to experimental vapor- 
liquid equilibrium data of pure ammonia. The resulting molecular model 
shows mean unsigned deviations to experiment of 0.7 % in saturated liq- 
uid density, 1.6 % in vapor pressure, and 2.7 % in enthalpy of vaporization 
over the whole temperature range from triple point to critical point. This 
new molecular model is used to predict thermophysical properties in the 
liquid, vapor and supercritical region, which are in excellent agreement 
with a high precision equation of state, that was optimized to 1147 ex- 
perimental data sets. Furthermore, it is also capable to predict the radial 
distribution functions properly, while no structural information is used in 
the optimization procedure. 
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24 1 Introduction 



25 Molecular modeling and simulation is a powerful tool for predicting thermo- 

26 physical properties, that is becoming more accesible due to the ever increasing 

27 computing power and the progress of methods and simulation tools. For real 

28 life applications in process engineering reliable predictions are needed for a wide 

29 variety of properties [Ij EJ . 

30 The central role for that task is played by the molecular model, that de- 

31 termines all of them. Therefore, a balanced modeling procedure, i.e. selection 

32 of model type and parameterization, is crucial. Unfortunately, thermophysi- 

33 cal properties usually depend on the model parameters in a highly non-linear 

34 fashion. So the development of new molecular models of technical quality is a 

35 time-consuming task. In this paper a procedure is proposed that uses informa- 

36 tion from ab initio quantum mechanical calculations to accelerate the modeling 

37 process. As an example, ammonia is regarded here. 

38 Ammonia is a well-known chemical intermediate, mostly used in fertilizer 

39 industries; another important application is its use as a refrigerant. Due to its 

40 simple symmetric structure and its strong intermolecular interactions it is also 

41 of high academic interest both experimentally and theoretically. 

42 Different approaches can be found in the literature to construct an inter- 

43 molecular potential for ammonia to be used in molecular simulation. Jorgensen 

44 and Ibrahim [?] as well as Hinchliffe et al. [5] used experimental bond distances 

45 and angles to place their interaction sites. Jorgensen and Ibrahim fitted a 12-6-3 

46 potential plus four partial charges to results from ab initio quantum mechanical 

47 calculations, they derived for 250 orientations of the ammonia dimer using the 

48 ST0-3G minimal basis set. To yield reasonable potential energies for liquid 

49 ammonia compared to experimental results, they had to scale their potential by 

50 a factor 1.26. 

51 Hinchliffe et al. used a combination of exponential repulsion terms, an at- 

52 tractive Morse potential, and four partial charges to construct the intermolecular 
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53 potential. The parameters were determined by fitting to a total of 61 points on 

54 the ammonia dimer energy surface at seven different orientations, which were 

55 calculated using the 6-31G* basis set. Hinchliffe et al. have pointed out, that 

56 the parameterization is ambiguous concerning the selection of dimer configu- 

57 rations and the used interaction potentials. Also the different models perform 

58 different well on various properties. 

59 In a later work Impey and Klein [6] reparameterized the molecular model 

60 by Hinchliffe et al. They switched to an "effective" pair potential using one 

61 Lennard- Jones (12-6) potential at the nitrogen nucleus site to describe the dis- 

62 persive and repulsive interactions. The parameters were optimized to the radial 

63 distribution function 5n-n of liquid ammonia measured by Narten [7]. 

64 Kristof et al. [S] used this model to predict vapor-liquid equilibrium prop- 

65 erties and found systematic deviations in both vapor pressure and saturated 

66 densities. So they decided to develop a completely new molecular model. Again 

67 they used experimental bond distances and angles to place the interaction sites. 

68 All further parameters of their model, i.e. the partial charges on all atoms and 

69 the parameters of the single Lennard- Jones (12-6) potential, were adjusted to 

70 vapor-liquid coexistence properties. With this model Kristof et al. reached 

71 a description of the vapor-liquid equilibrium (VLE) of ammonia of reasonable 

72 accuracy. 

73 For their simulations, Kristof et al. used the Gibbs ensemble Monte Carlo 

74 (GEMC) technique [3 [TO] with an extension to the NpH ensemble [HI [HI- This 

75 methods have some difficulties simulating strongly interacting fluids, yielding to 

76 relatively large statistical uncertainties. When applying our methods for the 

77 simulation of VLE, leading to much smaller statistical errors, we get results 

78 slightly outside the error bars of Kristof et al. Also systematic deviations to the 

79 experimental vapor-liquid properties are seen, especially in the vapor pressure 

80 and critical temperature, cf. Figures [T] to [31 

81 In the present work a new molecular model for ammonia is proposed. This 

82 model is based on the work of Kristof et al. and improved by including data on 
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83 geometry and electrostatics from ab initio quantum mechanical calculations. 

84 The paper is structured as follows: Initially, a procedure is proposed for the 

85 development of a preliminary molecular model. This model, called first model in 

86 the following, is then adjusted to experimental VLE data until a desired quality 

87 is reached. The resulting model, denoted as new model in the following, is 

88 used afterwards to predict thermal and caloric properties apart from the phase 

89 coexistence as well as structural properties. 



90 2 Selection of Model Type and Parameteriza- 

91 tion 

92 The modeling philosophy followed here is to keep the molecular model as simple 

93 as possible. Therefore, the molecule is assumed rigid and non-polarizable, i.e. 

94 a single state-independent set of parameters is used. Hydrogen atoms are not 

95 modeled explicitely, a united-atom approach is used. 

96 For both present models, a single Lennard- Jones potential was assumed to 

97 describe the dispersive and repulsive interactions. The electrostatic interactions 

98 as well as hydrogen bonding were modeled by a total of four partial charges. 

99 This modeling approach was found to be appropriate for other hydrogen bonding 

100 fluids like methanol T3] , ethanol P^l] , and formic acid [IS] and was also followed 

101 by Impey and Klein 6 and Kristof et al. 8^ for ammonia. 

102 Thus, the potential energy Uij between two ammonia molecules i and j is 

103 given by 

104 where a is the site index of charges on molecule i and b the site index of charges 

105 on molecule j, respectively. The site-site distances between molecules i and j 

106 are denoted by r.y for the single Lennard- Jones potential and rijab for the four 

107 partial charges, respectively, a and e are the Lennard-Jones size and energy 

108 parameters, while qia and qjb are the partial charges located at the sites a and 
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109 b on the molecules i and respectively. Finally, eg denotes the permittivity of 

no the vacuum. 

111 To keep the modeling procedure as independent as possible from the avail- 

112 ability of specific information, no experimental bond lengths or angles were used 

113 here in contrast to [H O [6l [8] . Instead, the nucleus positions were calculated by 

114 the means of quantum mechanics, where the software package GAMESS (US) 

115 |16j was used. A geometry optimization was performed on the Hartree-Fock, 

116 i.e. self-consistent field (SCF), level using the basis set 6-31G, which is a split- 

117 valence orbital basis set without polarizable terms. The nucleus positions from 
lis this ab initio calculation were directly used to specifiy the positions of the five 

119 interaction sites. At the nitrogen nucleus site and at each of the hydrogen nu- 

120 cleus sites a partial charge was placed. The Lennard-Jones site conincides with 

121 the nitrogen nucleus position, cf. Table [TJ 

122 To obtain the magnitude of the partial charges, another subsequent quantum 

123 mechanical calculation was performed. It was done on M0ller-Plesset 2 level 

124 using the polarizable basis set 6-311G(d,p) and the geometry from the previous 

125 step. By default, quantum mechanical calculations are performed on a single 

126 molecule of interest in vacuum. It is widely known, that the gas phase dipolar 

127 moments significantly differ from the dipole moment in the liquid state. As it 

128 was seen from former work [17[ I18j , molecular models yield better results on 

129 VLE properties, when a " liquid- like" dipolar moment is applied. Therefore, the 

130 single molecule was calculated within a dielectric cavity utilizing the COSMO 

131 (COnducter like Screening MOdel) method [12] to mimic the liquid state. The 

132 partial charges were chosen to yield the resulting dipole moment of 1.94 Debye, 

133 the parameters are given in Table [1] 

134 The first model combines this electrostatics with the Lennard-Jones parame- 

135 ters of Kristof et al [8J , so no additional experimental data was used. To achieve 

136 an optimized new model, the two Lennard-Jones parameters a and e were ad- 

137 justed to experimental saturated liquid density, vapor pressure, and enthalpy of 

138 vaporization using a Newton scheme as proposed by Stoll [20] . These properties 
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139 were chosen for the adjustment as they all represent major characteristics of 

140 the fluid region. Furthermore, they are relatively easy to be measured and are 

141 available for many components of technical interest. 

142 The applied optimization method has many similarities with the one pro- 

143 posed by lingerer et al. 21j and later on modified by Bourasseau et al. [22]. It 

144 relies on a least-square minimiztion of a weighted fitness function J- that quan- 

145 tifies the devitions of simulation results from a given molecular model compared 

146 to experimental data. The weighted fitness function writes as 

^ d ^ 

^=lYl 711 M,sUMo) - A,,exp)' , (2) 

147 wherein the n-dimensional vector Mq — {ttiq^i, ...,mo,„) is a short-cut notation 

148 for the set of n model parameters too,i, TOo.n to be optimized. The deviations 

149 of results from simulation Ai^sim to experimental data ^i.cxp are weighted with 

150 the expected simulation uncertainties SAi^sim- Equation ^ allows simultaneous 

151 adjustment of the model parameters to different thermophysical properties A 

152 (saturated liquid densities p', vapor pressures p'^, and enthalpies of vaporization 

153 A/iv at various temperatures in the present work). 

154 The unknown functional dependence of the property A on the model param- 

155 eters is approximated by a first order Taylor series developed in the vicinity of 

156 the initial parameter set Mq 

" QA, Sim 

Ai_sim(M„ow) = A^^s\ui{Mo) + 2_] Q^"^ ■ (Wncwj ^ moj) . (3) 

157 Therein, the partial derivatives of Ai with respect to each model parameter rrij, 

158 i.e. the sensitivities, are calculated from difference quotients 

dAj^si^ ^ Ai,sim("T-oa, ■■■,moj + Amj, ...,too^„) - >li,sim(wo4, ...,moj, ...,too^„) 
drrij Arrij 

(4) 

159 Assuming a sound choice of the model parameter variations Arrij, i.e. small 

160 enough to ensure linearity and large enough to yield differences in the simulation 

161 results significantly above the statistical uncertainties, this method allows a step- 

162 wise optimization of the molecular model by minimization of the fitness function 

163 J-. Experience shows that an optimized set of model parameters can be found 
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164 within a few steps when starting from a reasonable initial model. 

3 Results and Discussion 

166 3.1 Vapor-Liquid Equilibria 

167 VLE results for the new model are compared to data obtained from a reference 

168 quality equation of state (EOS) [53] in Figures [1] to [3] These figures also include 

169 the results, that we calculated using the first model and the model from Kristof 

170 et al. [5]. The present numerical simulation results together with experimental 

171 data |23j are given in Table [21 technical simulation details are given in the 

172 appendix. 

173 The reference EOS [23j used for adjustment and comparison here, was op- 

174 timized to 1147 experimental data sets. It is based on two older EOS from 

175 the late nineteen seventies [211 125] and also recommended by the NIST within 

176 their reference EOS database REFPROP [26^. The proposed uncertainties of 

177 the equation of state are 0.2 % in density, 2 % in heat capacity, and 2 % in the 

178 speed of sound, except in the critical region. The uncertainty in vapor pressure 

179 is 0.2 %. 

180 The model of Kristof et al. shows noticeable deviations from experimental 

181 data. The mean unsigned errors over the range of VLE are 1.9 % in saturated 

182 liquid density, 13 % in vapor pressure and 5.1 % in enthalpy of vaporization. 

183 Even without any further adjustment to experimental data a better description 

184 was found using the first model. The deviations between simulation results and 

185 reference EOS are 1.5 % in saturated liquid density, 10.4 % in vapor pressure 

186 and 5.1 % in enthalpy of vaporization. 

187 With the new model a significant improvement is achieved compared to the 

188 model from Kristof et al. The description of the experimental VLE is very 

189 good, the mean unsigned deviations in saturated liquid density, vapor pressure 

190 and enthalpy of vaporization are 0.7, 1.6, and 2.7 %, respectively. Only at low 

191 temperatures, in the range of ambient pressure, a slightly worse description of 
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192 the vapor pressure compared to the first model is observed. In Figure 2] the 

193 relative deviations of the new model, the model from Kristof et al., and the first 

194 model are shown in the whole range of the VLE starting from triple point to 

195 critical point. 

196 Mathews |27j gives experimental critical values of temperature, density and 

197 pressure for ammonia: Tc ==405.65 K, =13.8 mol/1, and pc =11.28 MPa. 

198 Following the procedure suggested by Lotfi et al. [28] the critical properties 

199 Tc =395.82 K, =14.0 mol/1, and pc =11.26 MPa for the model of Kristof et al. 

200 were calculated, where the critical temperature is underestimated by 2.4 %. For 

201 the first model Tc =403.99 K, pc =14.1 mol/1, andpc =11.67 MPa were obtained 

202 and for the new model Tc =402.21 K, pc =13.4 mol/1, and pc =10.52 MPa. The 

203 latter two give reasonable results for the critical temperature, while the new 

204 model underpredicts the critical pressure slightly. 

205 3.2 Homogeneous Region 

206 In many technical applications thermodynamic properties in the homogeneous 

207 fiuid region apart from the VLE are needed. Thus, the new molecular model 

208 was tested on its predictive capabilities in these states. 

209 Thermal and caloric properties were predicted with the new model in the 

210 homogenous liquid, vapor and supercritical fluid region. In total, 70 state points 

211 were regarded, covering a large range of states with temperatures up to 700 K 

212 and pressures up to 700 MPa. In Figure O relative deviations between simula- 

213 tion and reference EOS 23J in terms of density are shown. The deviations are 

214 typically below 3 % with the exception of the extended critical region, where a 

215 maximum deviation of 6.8 % is found. 

216 Figure[H]presents relative deviations in terms of enthalpy between simulation 

217 and reference EOS 23J. In this case deviations are very low for low pressures 

218 and high temperatures (below 1-2 %). Typical deviations in the other cases are 

219 below 5 %. 

220 These results confirm the modeling procedure. By adjustment to VLE data 
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221 only, quantitatively correct predictions in most of the technically important fluid 

222 region can be obtained. 

223 3.3 Second Virial Coefficient 

224 The virial expansion gives an equation of state for low density gases. For am- 

225 monia it is a good approximation for gaseous states below 0.1 MPa with a 

226 maximum error of 2.5 %. Starting from the nineteen thirties it was shown, that 

227 the virial coefficients can nowadays easily be derived from the intermolecular 

228 potential [221 [301 [SI] ■ The second virial coefhcient is related to the molecular 

229 model by [32] 

230 where (r^ , a;^, u;-,) is the interaction energy between two molecules i and j, 

231 cf. Equation ([T]). /cb denotes Boltzmann's constant and the () brackets indicate 

232 an average over the orientations u^i and u)j of the two molecules separated by 

233 the center of mass distance rij . 

234 The second virial coefficient was calculated here by evaluating Mayer's /- 

235 function at 363 radii from 2.4 to 8 A, averaging over 500^ random orientations 

236 at each radius. The random orientations were generated using a modified Monte 

237 Carlo scheme [331 [2 ■ A cut-off correction was applied for distances larger than 

238 8 A for the LJ potential [34 . The electrostatic interactions need no long-range 

239 correction as they vanish by angle averaging. 

240 Figure [7] shows the second virial coefficient predicted by the new model is 

241 shown in comparison to the reference EOS. An excellent agreement was found 

242 over the full temperature range with a maximum deviation of —4.3 % at 300 K. 

243 3.4 Structural Quantities 

244 Due to its scientific and technical importance, experimental data on the micro- 

245 scopic structure of liquid ammonia are available. Narten [^ and Ricci et al. |35| 

246 applied X-ray and neutron diffraction, respectively. The results from Ricci et 
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247 al. show a smoother gradient and are available for all three types of atom-atom 

248 pair correlations, namely nitrogen- nitrogen (N-N), nitrogen-hydrogen (N-H), 

249 and hydrogen- hydrogen (H-H), thus they were used for comparison here. In 

250 Figure [8l these experimental radial distribution functions for liquid ammonia 

251 at 273.15 K and 0.483 MPa are compared to present predictive simulation data 

252 based on the new model. 

253 It is found that these structural properties are in very good agreement, 

254 although no adjustment was done with regard to structural properties. The 

255 atom-atom distance of the first three layers is predicted correctly, while only 

256 minor overshootings in the first peak are found. Please note, that the first peak 

257 of experiment in 5n-h and .gn-H show intramolecular pair correlations, which 

258 are not calculated in the simulation. 

259 In the experimental radial distribution function .gN-H the hydrogen bonding 

260 of ammonia can be seen at 2-2.5 A. Due to the simplified approximation by 

261 off-centric partial charges, the molecular model is not capable to describe this 

262 effect completely. But even with this simple model a small shoulder at 2.5 A is 

263 obtained. 

264 4 Conclusion 

265 A new molecular model is proposed for ammonia. This model was developed 

266 using a new modeling procedure, which speeds up the modeling process and can 

267 be applied on arbitrary molecules. The interaction sites are located according to 

268 atom positions resulting from ab initio quantum mechanical calculations. Also 

269 the electrostatic interactions, here in form of partial charges, are parameterized 

270 according to high-level ab initio quantum mechanical results. The latter are 

271 obtained by calculations within a dielectric continuum to mimic the (stronger) 

272 interactions in the liquid phase. The partial charges for the present ammonia 

273 model are specified to yield the same dipole moment as quantum mechanics. The 

274 Lennard-Jones parameters were adjusted to VLE data, namely vapor pressure. 
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275 saturated liquid density, and enthalpy of vaporization. 

276 A description of the VLE of ammonia was reached within relative deviations 

277 of a few percents. Next to this, covering a large region of states, a good pre- 

278 diction of both thermal and caloric properties apart from the VLE was found 

279 compared to a reference EOS [23]. 

280 Predicted structural quantities, i.e. radial distribution functions in the liquid 

281 state, are in very good agreement to experimental neutron diffraction data. 

282 This shows, that molecular models adjusted to macroscopic thermodynamic 

283 properties also give reasonable results on microscopic properties. Note that this 

284 is not true vice versa in most cases. With the present model a similar quality 

285 in describing the atomic radial distribution functions as Impey and Klein [S] is 

286 gained, while the macroscopic properties like vapor pressure differ considerably. 

287 So the latter can be seen as the more demanding criteria. 
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295 6 Appendix 

296 The Grand Equilibrium method ^36j was used to calculate VLE data at eight 

297 temperatures from 240 to 395 K during the optimization process. At each 

298 temperature for the liquid, molecular dynamics simulations were performed in 

299 the NpT ensemble using isokinetic velocity scaling 34j and Anderson's barostat 

300 [27]. There, the number of molecules is 864 and the time step was 0.58 fs except 
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301 for the lowest temperature, where 1372 molecules and a time step of 0.44 fs 

302 were used. The initial configuration was a face centered cubic lattice, the fluid 

303 was equilibrated over 120 000 time steps with the first 20 000 time steps in the 

304 canonical (NVT) ensemble. The production run went over 300 000 time steps 

305 (400 000 for 240 K) with a membrane mass of 10^ kg/m*. Widom's insertion 

306 method j38j was used to calculate the chemical potential by inserting up to 4 000 

307 test molecules every production time step. 

308 At the lowest two temperatures additional Monte Carlo simulations were 

309 performed in the NpT ensemble for the liquid. There, the chemical potential of 

310 liquid ammonia was calculated by the gradual insertion method [40J . The num- 

311 ber of molecules was 500. Starting from a face centered cubic lattice, 15 000 

312 Monte Carlo cycles were performed for equilibration and 50 000 for production, 

313 each cycle containing 500 displacement moves, 500 rotation moves, and 1 volume 

314 move. Every 50 cycles 5000 fluctuating state change moves, 5000 fluctuating par- 

315 tide translation/rotation moves, and 25000 biased particle translation/rotation 

316 moves were performed, to measure the chemical potential. These computation- 

317 ally demanding simulations yield the chemical potential in the dense and strong 

318 interacting liquid with high accuracy, leading to small uncertainties in the VLE. 

319 For the corresponding vapor, Monte Carlo simulations in the psendo-fiVT 

320 ensemble were performed. The simulation volume was adjusted to lead to an 

321 average number of 500 molecules in the vapor phase. After 1 000 initial NVT 

322 Monte Carlo cycles, starting from a face centered cubic lattice, 10 000 equi- 

323 libration cycles in the pseudo-/xyr ensemble were performed. The length of 

324 the production run was 50 000 cycles. One cycle is defined here to be a num- 

325 ber of attempts to displace and rotate molecules equal to the actual number of 

326 molecules plus three insertion and three deletion attempts. 

327 The cut-off radius was set to 17.5 A throughout and a center of mass cut-off 

328 scheme was employed. Lennard- Jones long-range interactions beyond the cut-off 

329 radius were corrected as proposed in [34 . Electrostatic interactions were ap- 

330 proximated by a resulting molecular dipole and corrected using the reaction field 
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331 method |34| . Statistical uncertainties in the simulated values were estimated by 

332 a block averaging method [H] . 

333 For the simulations in the homogeneous region, molecular dynamics sim- 

334 ulations were performed with the same technical parameters as used for the 

335 saturated liquid runs. 

336 For the radial distribution functions a molecular dynamics simulation was 

337 performed with 500 molecules. Intermolecular site-site distances were divided 

338 in 200 slabs from to 13.5 A and summed up for 50 000 time steps. 
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Table 1; Parameters of the new ammonia model. The electronic charge is 
e = 1.6021 • 10-19 C. 



Interaction 


X 


y 


z 


a 




9 


Site 


A 


A 


A 


A 


K 


e 


N 


0.0 


0.0 


0.0757 


3.376 


182.9 


-0.9993 


H(l) 


6.9347 


6.0 


-0.3164 






6.3331 


H(2) 


-0.4673 


6.8095 


-0.3164 






6.3331 


H(3) 


-0.4673 


-0.8095 


-0.3164 






6.3331 
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Table 2: Vapor- liquid equilibria of ammonia: simulation results using the new 
model (sim) compared to data from a reference quality equation of state [25] 
(eos) for vapor pressure, saturated densities and enthalpy of vaporization. The 
number in parentheses indicates the statistical uncertainty in the last digit. 



T 


Psim 


Pcos 


Pshn 


Pcos 


p"- 

r Sim 


Pcos 






K 


MPa 


MPa 


mol/l 


mol/l 


mol/l 


mol/l 


k J / mol 


k J / mol 


240 


0.12(1) 


0.102 


40.26(1) 


40.032 


0.066(5) 


0.0527 


24.11(1) 


23.31 


280 


0.60(2) 


0.551 


36.98(2) 


36.939 


0.280(8) 


0.257 


21.56(1) 


21.07 


315 


1.65(4) 


1.637 


33.76(3) 


33.848 


0.74 (1) 


0.744 


18.96(2) 


18.57 


345 


3.37(4) 


3.457 


30.45(4) 


30.688 


1.55 (1) 


1.624 


16.19(3) 


15.79 


363 


5.22(5) 


5.101 


28.17(6) 


28.368 


2.56 (2) 


2.544 


13.93(5) 


13.65 


375 


6.37(6) 


6.485 


26.18(7) 


26.502 


3.17 (3) 


3.459 


12.48(6) 


11.89 


385 


7.88(5) 


7.845 


24.05(9) 


24.608 


4.27 (5) 


4.554 


10.49(9) 


10.08 


395 


9.54(7) 


9.422 


20.9 (1) 


22.090 


5.66 (9) 


6.272 


8.1 (1) 


7.66 
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Abstract 

An optimized molecular model for ammonia, which is based on a previ- 
ous work of Kristof et al., Mol. Phys. 97 (1999) 1129-1137, is presented. 
Improvements are achieved by including data on geometry and electro- 
statics from ab initio quantum mechanical calculations in a first model. 
Afterwards the parameters of the Lennard- Jones potential, modeling dis- 
persive and repulsive interactions, are optimized to experimental vapor- 
liquid equilibrium data of pure ammonia. The resulting molecular model 
shows mean unsigned deviations to experiment of 0.7 % in saturated liq- 
uid density, 1.6 % in vapor pressure, and 2.7 % in enthalpy of vaporization 
over the whole temperature range from triple point to critical point. This 
new molecular model is used to predict thermophysical properties in the 
liquid, vapor and supercritical region, which are in excellent agreement 
with a high precision equation of state, that was optimized to 1147 ex- 
perimental data sets. Furthermore, it is also capable to predict the radial 
distribution functions properly, while no structural information is used in 
the optimization procedure. 
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24 1 Introduction 



25 Molecular modeling and simulation is a powerful tool for predicting thermo- 

26 physical properties, that is becoming more accesible due to the ever increasing 

27 computing power and the progress of methods and simulation tools. For real 

28 life applications in process engineering reliable predictions are needed for a wide 

29 variety of properties [Ij EJ . 

30 The central role for that task is played by the molecular model, that de- 

31 termines all of them. Therefore, a balanced modeling procedure, i.e. selection 

32 of model type and parameterization, is crucial. Unfortunately, thermophysi- 

33 cal properties usually depend on the model parameters in a highly non-linear 

34 fashion. So the development of new molecular models of technical quality is a 

35 time-consuming task. In this paper a procedure is proposed that uses informa- 

36 tion from ab initio quantum mechanical calculations to accelerate the modeling 

37 process. As an example, ammonia is regarded here. 

38 Ammonia is a well-known chemical intermediate, mostly used in fertilizer 

39 industries; another important application is its use as a refrigerant. Due to its 

40 simple symmetric structure and its strong intermolecular interactions it is also 

41 of high academic interest both experimentally and theoretically. 

42 Different approaches can be found in the literature to construct an inter- 

43 molecular potential for ammonia to be used in molecular simulation. Jorgensen 

44 and Ibrahim [?] as well as Hinchliffe et al. [5] used experimental bond distances 

45 and angles to place their interaction sites. Jorgensen and Ibrahim fitted a 12-6-3 

46 potential plus four partial charges to results from ab initio quantum mechanical 

47 calculations, they derived for 250 orientations of the ammonia dimer using the 

48 ST0-3G minimal basis set. To yield reasonable potential energies for liquid 

49 ammonia compared to experimental results, they had to scale their potential by 

50 a factor 1.26. 

51 Hinchliffe et al. used a combination of exponential repulsion terms, an at- 

52 tractive Morse potential, and four partial charges to construct the intermolecular 
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53 potential. The parameters were determined by fitting to a total of 61 points on 

54 the ammonia dimer energy surface at seven different orientations, which were 

55 calculated using the 6-31G* basis set. Hinchliffe et al. have pointed out, that 

56 the parameterization is ambiguous concerning the selection of dimer configu- 

57 rations and the used interaction potentials. Also the different models perform 

58 different well on various properties. 

59 In a later work Impey and Klein [6] reparameterized the molecular model 

60 by Hinchliffe et al. They switched to an "effective" pair potential using one 

61 Lennard- Jones (12-6) potential at the nitrogen nucleus site to describe the dis- 

62 persive and repulsive interactions. The parameters were optimized to the radial 

63 distribution function 5n-n of liquid ammonia measured by Narten [7]. 

64 Kristof et al. [S] used this model to predict vapor-liquid equilibrium prop- 

65 erties and found systematic deviations in both vapor pressure and saturated 

66 densities. So they decided to develop a completely new molecular model. Again 

67 they used experimental bond distances and angles to place the interaction sites. 

68 All further parameters of their model, i.e. the partial charges on all atoms and 

69 the parameters of the single Lennard- Jones (12-6) potential, were adjusted to 

70 vapor-liquid coexistence properties. With this model Kristof et al. reached 

71 a description of the vapor-liquid equilibrium (VLE) of ammonia of reasonable 

72 accuracy. 

73 For their simulations, Kristof et al. used the Gibbs ensemble Monte Carlo 

74 (GEMC) technique [3 [TO] with an extension to the NpH ensemble [HI [HI- This 

75 methods have some difficulties simulating strongly interacting fluids, yielding to 

76 relatively large statistical uncertainties. When applying our methods for the 

77 simulation of VLE, leading to much smaller statistical errors, we get results 

78 slightly outside the error bars of Kristof et al. Also systematic deviations to the 

79 experimental vapor-liquid properties are seen, especially in the vapor pressure 

80 and critical temperature, cf. Figures [T] to [31 

81 In the present work a new molecular model for ammonia is proposed. This 

82 model is based on the work of Kristof et al. and improved by including data on 
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83 geometry and electrostatics from ab initio quantum mechanical calculations. 

84 The paper is structured as follows: Initially, a procedure is proposed for the 

85 development of a preliminary molecular model. This model, called first model in 

86 the following, is then adjusted to experimental VLE data until a desired quality 

87 is reached. The resulting model, denoted as new model in the following, is 

88 used afterwards to predict thermal and caloric properties apart from the phase 

89 coexistence as well as structural properties. 



90 2 Selection of Model Type and Parameteriza- 

91 tion 

92 The modeling philosophy followed here is to keep the molecular model as simple 

93 as possible. Therefore, the molecule is assumed rigid and non-polarizable, i.e. 

94 a single state-independent set of parameters is used. Hydrogen atoms are not 

95 modeled explicitely, a united-atom approach is used. 

96 For both present models, a single Lennard- Jones potential was assumed to 

97 describe the dispersive and repulsive interactions. The electrostatic interactions 

98 as well as hydrogen bonding were modeled by a total of four partial charges. 

99 This modeling approach was found to be appropriate for other hydrogen bonding 

100 fluids like methanol T3] , ethanol P^l] , and formic acid [IS] and was also followed 

101 by Impey and Klein 6 and Kristof et al. 8^ for ammonia. 

102 Thus, the potential energy Uij between two ammonia molecules i and j is 

103 given by 

104 where a is the site index of charges on molecule i and b the site index of charges 

105 on molecule j, respectively. The site-site distances between molecules i and j 

106 are denoted by r.y for the single Lennard- Jones potential and rijab for the four 

107 partial charges, respectively, a and e are the Lennard-Jones size and energy 

108 parameters, while qia and qjb are the partial charges located at the sites a and 
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109 b on the molecules i and respectively. Finally, eg denotes the permittivity of 

no the vacuum. 

111 To keep the modeling procedure as independent as possible from the avail- 

112 ability of specific information, no experimental bond lengths or angles were used 

113 here in contrast to [H O [6l [8] . Instead, the nucleus positions were calculated by 

114 the means of quantum mechanics, where the software package GAMESS (US) 

115 |16j was used. A geometry optimization was performed on the Hartree-Fock, 

116 i.e. self-consistent field (SCF), level using the basis set 6-31G, which is a split- 

117 valence orbital basis set without polarizable terms. The nucleus positions from 
lis this ab initio calculation were directly used to specifiy the positions of the five 

119 interaction sites. At the nitrogen nucleus site and at each of the hydrogen nu- 

120 cleus sites a partial charge was placed. The Lennard-Jones site conincides with 

121 the nitrogen nucleus position, cf. Table [TJ 

122 To obtain the magnitude of the partial charges, another subsequent quantum 

123 mechanical calculation was performed. It was done on M0ller-Plesset 2 level 

124 using the polarizable basis set 6-311G(d,p) and the geometry from the previous 

125 step. By default, quantum mechanical calculations are performed on a single 

126 molecule of interest in vacuum. It is widely known, that the gas phase dipolar 

127 moments significantly differ from the dipole moment in the liquid state. As it 

128 was seen from former work [17[ I18j , molecular models yield better results on 

129 VLE properties, when a " liquid- like" dipolar moment is applied. Therefore, the 

130 single molecule was calculated within a dielectric cavity utilizing the COSMO 

131 (COnducter like Screening MOdel) method [12] to mimic the liquid state. The 

132 partial charges were chosen to yield the resulting dipole moment of 1.94 Debye, 

133 the parameters are given in Table [1] 

134 The first model combines this electrostatics with the Lennard-Jones parame- 

135 ters of Kristof et al [8J , so no additional experimental data was used. To achieve 

136 an optimized new model, the two Lennard-Jones parameters a and e were ad- 

137 justed to experimental saturated liquid density, vapor pressure, and enthalpy of 

138 vaporization using a Newton scheme as proposed by Stoll [20] . These properties 
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139 were chosen for the adjustment as they all represent major characteristics of 

140 the fluid region. Furthermore, they are relatively easy to be measured and are 

141 available for many components of technical interest. 

142 The applied optimization method has many similarities with the one pro- 

143 posed by lingerer et al. 21j and later on modified by Bourasseau et al. [22]. It 

144 relies on a least-square minimiztion of a weighted fitness function J- that quan- 

145 tifies the devitions of simulation results from a given molecular model compared 

146 to experimental data. The weighted fitness function writes as 

^ d ^ 

^=lYl 711 M,sUMo) - A,,exp)' , (2) 

147 wherein the n-dimensional vector Mq — {ttiq^i, ...,mo,„) is a short-cut notation 

148 for the set of n model parameters too,i, TOo.n to be optimized. The deviations 

149 of results from simulation Ai^sim to experimental data ^i.cxp are weighted with 

150 the expected simulation uncertainties SAi^sim- Equation ^ allows simultaneous 

151 adjustment of the model parameters to different thermophysical properties A 

152 (saturated liquid densities p', vapor pressures p'^, and enthalpies of vaporization 

153 A/iv at various temperatures in the present work). 

154 The unknown functional dependence of the property A on the model param- 

155 eters is approximated by a first order Taylor series developed in the vicinity of 

156 the initial parameter set Mq 

" QA, Sim 

Ai_sim(M„ow) = A^^s\ui{Mo) + 2_] Q^"^ ■ (Wncwj ^ moj) . (3) 

157 Therein, the partial derivatives of Ai with respect to each model parameter rrij, 

158 i.e. the sensitivities, are calculated from difference quotients 

dAj^si^ ^ Ai,sim("T-oa, ■■■,moj + Amj, ...,too^„) - >li,sim(wo4, ...,moj, ...,too^„) 
drrij Arrij 

(4) 

159 Assuming a sound choice of the model parameter variations Arrij, i.e. small 

160 enough to ensure linearity and large enough to yield differences in the simulation 

161 results significantly above the statistical uncertainties, this method allows a step- 

162 wise optimization of the molecular model by minimization of the fitness function 

163 J-. Experience shows that an optimized set of model parameters can be found 
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164 within a few steps when starting from a reasonable initial model. 

3 Results and Discussion 

166 3.1 Vapor-Liquid Equilibria 

167 VLE results for the new model are compared to data obtained from a reference 

168 quality equation of state (EOS) [53] in Figures [1] to [3] These figures also include 

169 the results, that we calculated using the first model and the model from Kristof 

170 et al. [5]. The present numerical simulation results together with experimental 

171 data |23j are given in Table [21 technical simulation details are given in the 

172 appendix. 

173 The reference EOS [23j used for adjustment and comparison here, was op- 

174 timized to 1147 experimental data sets. It is based on two older EOS from 

175 the late nineteen seventies [211 125] and also recommended by the NIST within 

176 their reference EOS database REFPROP [26^. The proposed uncertainties of 

177 the equation of state are 0.2 % in density, 2 % in heat capacity, and 2 % in the 

178 speed of sound, except in the critical region. The uncertainty in vapor pressure 

179 is 0.2 %. 

180 The model of Kristof et al. shows noticeable deviations from experimental 

181 data. The mean unsigned errors over the range of VLE are 1.9 % in saturated 

182 liquid density, 13 % in vapor pressure and 5.1 % in enthalpy of vaporization. 

183 Even without any further adjustment to experimental data a better description 

184 was found using the first model. The deviations between simulation results and 

185 reference EOS are 1.5 % in saturated liquid density, 10.4 % in vapor pressure 

186 and 5.1 % in enthalpy of vaporization. 

187 With the new model a significant improvement is achieved compared to the 

188 model from Kristof et al. The description of the experimental VLE is very 

189 good, the mean unsigned deviations in saturated liquid density, vapor pressure 

190 and enthalpy of vaporization are 0.7, 1.6, and 2.7 %, respectively. Only at low 

191 temperatures, in the range of ambient pressure, a slight deterioration compared 



7 



192 to the initial model is observed. In Figure [3] the relative deviations of the new 

193 model, the model from Kristof et al., and the first model are shown in the whole 

194 range of the VLE starting from triple point to critical point. 

195 Mathews |27j gives experimental critical values of temperature, density and 

196 pressure for ammonia: Tc =405.65 K, =13.8 mol/1, and pc =11.28 MPa. 

197 Following the procedure suggested by Lotfi et al. ^8] the critical properties 

198 Tc =395.82 K, =14.0 mol/1, and pc =11.26 MPa for the model of Kristof et al. 

199 were calculated, where the critical temperature is underestimated by 2.4 %. For 

200 the first model Tc =403.99 K, pc =14.1 mol/1, and pc =11.67 MPa were obtained 

201 and for the new model Tc =402.21 K, p^ =13.4 mol/1, and p^ =10.52 MPa. The 

202 latter two give reasonable results for the critical temperature, while the new 

203 model underpredicts the critical pressure slightly. 

204 3.2 Homogeneous Region 

205 In many technical applications thermodynamic properties in the homogeneous 

206 fluid region apart from the VLE are needed. Thus, the new molecular model 

207 was tested on its predictive capabilities in these states. 

208 Thermal and caloric properties were predicted with the new model in the 

209 homogenous liquid, vapor and supercritical fluid region. In total, 70 state points 

210 were regarded, covering a large range of states with temperatures up to 700 K 

211 and pressures up to 700 MPa. In Figure relative deviations between simula- 

212 tion and reference EOS 23J in terms of density are shown. The deviations are 

213 typically below 3 % with the exception of the extended critical region, where a 

214 maximum deviation of 6.8 % is found. 

215 Figure[H]presents relative deviations in terms of enthalpy between simulation 

216 and reference EOS 23J. In this case deviations are very low for low pressures 

217 and high temperatures (below 1-2 %). Typical deviations in the other cases are 

218 below 5 %. 

219 These results confirm the modeling procedure. By adjustment to VLE data 

220 only, quantitatively correct predictions in most of the technically important fluid 
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221 region can be obtained. 

222 3.3 Second Virial Coefficient 

223 The virial expansion gives an equation of state for low density gases. For am- 

224 monia it is a good approximation for gaseous states below 0.1 MPa with a 

225 maximum error of 2.5 %. Starting from the nineteen thirties it was shown, that 

226 the virial coefficients can nowadays easily be derived from the intermolecular 

227 potential [IHl 1301 131] ■ The second virial coefficient is related to the molecular 

228 model by [32] 

229 where Uij{rij,u}i,ujj) is the interaction energy between two molecules i and j, 

230 cf. Equation ([T]). denotes Boltzmann's constant and the () brackets indicate 

231 an average over the orientations <jJi and u}j of the two molecules separated by 

232 the center of mass distance r^j . 

233 The second virial coefficient was calculated here by evaluating Mayer's /- 

234 function at 363 radii from 2.4 to 8 A, averaging over 500^ random orientations 

235 at each radius. The random orientations were generated using a modified Monte 

236 Carlo scheme [331 12] ■ A cut-off correction was applied for distances larger than 

237 8 A for the LJ potential [34 . The electrostatic interactions need no long-range 

238 correction as they vanish by angle averaging. 

239 Figure [7] shows the second virial coefficient predicted by the new model is 

240 shown in comparison to the reference EOS. An excellent agreement was found 

241 over the full temperature range with a maximum deviation of —4.3 % at 300 K. 

242 3.4 Structural Quantities 

243 Due to its scientific and technical importance, experimental data on the micro- 

244 scopic structure of liquid ammonia are available. Narten [7] and Ricci et al. [35] 

245 applied X-ray and neutron diffraction, respectively. The results from Ricci et 

246 al. show a smoother gradient and are available for all three types of atom-atom 
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247 pair correlations, namely nitrogen- nitrogen (N-N), nitrogen-hydrogen (N-H), 

248 and hydrogen- hydrogen (H-H), thus they were used for comparison here. In 

249 Figure [51 these experimental radial distribution functions for liquid ammonia 

250 at 273.15 K and 0.483 MPa are compared to present predictive simulation data 

251 based on the new model. 

262 It is found that these structural properties are in very good agreement, 

253 although no adjustment was done with regard to structural properties. The 

254 atom-atom distance of the first three layers is predicted correctly, while only 

255 minor overshootings in the first peak are found. Please note, that the first peak 

256 of experiment in ^n-h and ^h-h show intramolecular pair correlations, which 

257 are not calculated in the simulation. 

258 In the experimental radial distribution function .gN-H the hydrogen bonding 

259 of ammonia can be seen at 2-2.5 A. Due to the simplified approximation by 

260 off-centric partial charges, the molecular model is not capable to describe this 

261 effect completely. But even with this simple model a small shoulder at 2.5 A is 

262 obtained. 

263 4 Conclusion 

264 A new molecular model is proposed for ammonia. This model was developed 

265 using a new modeling procedure, which speeds up the modeling process and can 

266 be applied on arbitrary molecules. The interaction sites are located according to 

267 atom positions resulting from ab initio quantum mechanical calculations. Also 

268 the electrostatic interactions, here in form of partial charges, are parameterized 

269 according to high-level ab initio quantum mechanical results. The latter are 

270 obtained by calculations within a dielectric continuum to mimic the (stronger) 

271 interactions in the liquid phase. The partial charges for the present ammonia 

272 model are specified to yield the same dipole moment as quantum mechanics. The 

273 Lennard-Jones parameters were adjusted to VLB data, namely vapor pressure, 

274 saturated liquid density, and enthalpy of vaporization. 
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275 A description of the VLE of ammonia was reached within relative deviations 

276 of a few percents. Next to this, covering a large region of states, a good pre- 

277 diction of both thermal and caloric properties apart from the VLE was found 

278 compared to a reference EOS [23]. 

279 Predicted structural quantities, i.e. radial distribution functions in the liquid 

280 state, are in very good agreement to experimental neutron diffraction data. 

281 This shows, that molecular models adjusted to macroscopic thermodynamic 

282 properties also give reasonable results on microscopic properties. Note that this 

283 is not true vice versa in most cases. With the present model a similar quality 

284 in describing the atomic radial distribution functions as Impey and Klein [B] is 

285 gained, while the macroscopic properties like vapor pressure differ considerably. 

286 So the latter can be seen as the more demanding criteria. 
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294 6 Appendix 

295 The Grand Equilibrium method 36] was used to calculate VLE data at eight 

296 temperatures from 240 to 395 K during the optimization process. At each 

297 temperature for the liquid, molecular dynamics simulations were performed in 

298 the NpT ensemble using isokinetic velocity scaling [34J and Anderson's barostat 

299 [37]. There, the number of molecules is 864 and the time step was 0.58 fs except 

300 for the lowest temperature, where 1372 molecules and a time step of 0.44 fs 
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301 were used. The initial configuration was a face centered cubic lattice, the fluid 

302 was equilibrated over 120 000 time steps with the first 20 000 time steps in the 

303 canonical (NVT) ensemble. The production run went over 300 000 time steps 

304 (400 000 for 240 K) with a membrane mass of 10^ kg/m*. Widom's insertion 

305 method [38] was used to calculate the chemical potential by inserting up to 4 000 

306 test molecules every production time step. 

307 At the lowest two temperatures additional Monte Carlo simulations were 

308 performed in the NpT ensemble for the liquid. There, the chemical potential of 

309 liquid ammonia was calculated by the gradual insertion method [40J . The num- 

310 ber of molecules was 500. Starting from a face centered cubic lattice, 15 000 

311 Monte Carlo cycles were performed for equilibration and 50 000 for production, 

312 each cycle containing 500 displacement moves, 500 rotation moves, and 1 volume 

313 move. Every 50 cycles 5000 fluctuating state change moves, 5000 fluctuating par- 

314 tide translation/rotation moves, and 25000 biased particle translation/rotation 

315 moves were performed, to measure the chemical potential. These computation- 

316 ally demanding simulations yield the chemical potential in the dense and strong 

317 interacting liquid with high accuracy, leading to small uncertainties in the VLE. 

318 For the corresponding vapor, Monte Carlo simulations in the psendo-iiVT 

319 ensemble were performed. The simulation volume was adjusted to lead to an 

320 average number of 500 molecules in the vapor phase. After 1 000 initial NVT 

321 Monte Carlo cycles, starting from a face centered cubic lattice, 10 000 equi- 

322 libration cycles in the paenAo-^VT ensemble were performed. The length of 

323 the production run was 50 000 cycles. One cycle is defined here to be a num- 

324 ber of attempts to displace and rotate molecules equal to the actual number of 

325 molecules plus three insertion and three deletion attempts. 

326 The cut-off radius was set to 17.5 A throughout and a center of mass cut-off 

327 scheme was employed. Lennard- Jones long-range interactions beyond the cut-off 

328 radius were corrected as proposed in [34 . Electrostatic interactions were ap- 

329 proximated by a resulting molecular dipole and corrected using the reaction field 

330 method [34] . Statistical uncertainties in the simulated values were estimated by 
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331 a block averaging method [H] . 

332 For the simulations in the homogeneous region, molecular dynamics sim- 

333 ulations were performed with the same technical parameters as used for the 

334 saturated liquid runs. 

335 For the radial distribution functions a molecular dynamics simulation was 

336 performed with 500 molecules. Intermolecular site-site distances were divided 

337 in 200 slabs from to 13.5 A and summed up for 50 000 time steps. 



13 



338 References 



339 [1] P. Ungerer, V. Lachet, and B. Tavitian. Applications of molecular simu- 

340 lation in oil and gas production and processing. Oil Gas Sci. TechnoL, 61 

341 (2006) 387-403. 

342 [2] B. Eckl, J. Vrabec, and H. Hasse. On the application of force fields for 

343 predicting a wide variety of properties: Ethylene oxide as an example. 

344 Fluid Phase Equilib., submitted (2007). 

345 [3] B. Eckl, Y.-L. Huang, J. Vrabec, and H. Hasse. Vapor pressure of R227ea 

346 + Ethanol at 343.17 K by molecular simulation. Fluid Phase Equilib., 260 

347 (2007) 177-182. 

348 [4] W.L. Jorgensen and M. Ibrahim. Structure and properties of liquid ammo- 

349 nia. J. Amer. Chem. Soc. 102 (1980) 3309-3315. 

350 [5] A. Hinchhffe, D.G. Bounds, M.L. Klein, I.R. McDonald, and R. Righini. 

351 Intermolecular potentials for ammonia based on SCF-MO calculations. J. 

352 Chem. Phys. 74 (1981) 1211-1216. 

353 [6] R.W. Impey and M.L. Klein. A simple intermolecular potential for liquid 

354 ammonia. Chem. Phys. Lett. 104 (1984) 579-582. 

355 [7] A.H. Narten. Liquid- ammonia - molecular correlation-functions from x- 

356 ray-diffraction. J. Chem. Phys. 66 (1977) 3117-3120. 

357 [8] T. Kristof, J. Vorholz, J. Liszi, B. Rumpf, and G. Maurer. A simple effective 

358 pair potential for the molecular simulation of the thermodynamic properties 

359 of ammonia. Mol. Phys. 97 (1999) 1129-1137. 

360 [9] A.Z. Panagiotopoulos. Direct determination of phase coexistence properties 

361 of fluids by Monte-Carlo simulation in a new ensemble. Mol. Phys., 61 

362 (1987) 813-826. 



14 



363 [10] A.Z. Panagiotopoulos, N. Quirke, M. Stapleton, and D.J. Tildesley. Phase- 

364 equilibria by simulation in the Gibbs ensemble - alternative derivation, 

365 generalization and application to mixture and membrane equilibria. Mol. 

366 Phys. 63 (1988) 527-545. 

367 [11] T. Kristof, J. Liszi. Alternative implementations of the Gibbs ensemble 

368 Monte Carlo calculation. Chem. Phys. Lett., 261 (1996) 620-624. 

369 [12] T. Kristof, J. Liszi. Application of a new Gibbs ensemble Monte Carlo 

370 method to site-site interaction model fluids. Mol. Phys., 90 (1997) 1031- 

371 1034. 

372 [13] T. Schnabcl, M. Cortada, J. Vrabec, S. Lago, and H. Hasse. Molecular 

373 model for formic acid adjusted to vapor-liquid equilibria. Chem. Phys. 

374 Lett., 435 (2007) 268-272. 

375 [14] T. Schnabel, J. Vrabec, and H. Hasse. Henry's law constants of methane, 

376 nitrogen, oxygen and carbon dioxide in ethanol from 273 to 498 K: predic- 

377 tion from molecular simulation. Fluid Phase Equilih., 233 (2005) 134-143. 

378 [15] T. Schnabel, A. Srivastava, J. Vrabec, and H. Hasse. Hydrogen Bonding 

379 of Methanol in Supercritical C02: Comparison between IH-NMR Spec- 

380 troscopic Data and Molecular Simulation Results, J. Phys. Chem. B, 111 

381 (2007) 9871-9878. 

382 [16] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, 

383 J.H. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, Kiet. General atomic 

384 and molecular electronic structure system. J. Comput. Chem., 14 (1993) 

385 1347-1363. 

386 [17] J. Vrabec, J. StoU, and H. Hasse. A set of molecular models for symmetric 

387 quadrupolar fluids. J. Phys. Chem. B, 105 (2001) 12126-12133. 



15 



388 [18] J. StoU, J. Vrabec, and H. Hasse. A set of molecular models for car- 

389 bon monoxide and halogenated hydrocarbons. J. Chem. Phys., 119 (2003) 

390 11396-11407. 

391 [19] K. Baldridge and A. Klamt. First principles implementation of solvent 

392 effects without outlying charge error. J. Chem. Phys., 106 (1997) 6622- 

393 6633. 

394 [20] J. StoU. Molecular Models for the Prediction of Thermophysical Properties 

395 of Pure Fluids and Mixtures. Fortschritt-Berichte VDI, Reihe 3, 836, VDI- 

396 Verlag, Diisseldorf, 2005. 

397 [21] P. Ungerer, A. Boutin, and A.H. Fuchs. Direct calculation of bubble points 

398 by Monte Carlo simulation. Mol. Phys., 97 (1999) 523-539. 

399 [22] E. Bourasseau, M. Haboudou, A. Boutin, A.H. Fuchs, and P. Ungerer. New 

400 optimization method for intermolecular potentials: Optimization of a new 

401 anisotropic united atoms potential for olefins: Prediction of equilibrium 

402 properties. J. Chem. Phys., 118 (2003) 3020-3034. 

403 [23] R. Tillner-Roth, F. Harms- Watzenberg, and H. D. Baehr. Eine neue Funda- 

404 mentalgleichung fuer Ammoniak. DKV-Tagungsbericht, 20 (1993) 167-181. 

405 [24] L. Haar. Thermodynamic properties of ammonia. J. Phys. Chem. Ref. 

406 Data, 7 (1978) 635-792. 

407 [25] J. Ahrendts and H.D. Baehr. Die thermodynamischen Eigenschaften von 

408 Ammoniak. VDI-Forschungsheft Nr. 596, VDI- Verlag, Diisseldorf, 1979. 

409 [26] E.W. Lemmon, M.L. Huber, and M.O. McLinden. REFPROP, NIST Stan- 
no dard Reference Database 23, Version 8.0. Physical and Chemical Properties 

411 Division, National Institute of Standards and Technology, Boulder, 200x. 

412 [27] J.F. Mathews. The critical constants of inorganic substances. Chem. Rev., 

413 72 (1972) 71-100. 



16 



414 [28] A. Lotfi, J. Vrabec, and J. Fischer, Vapour liquid equilibria of the Lennard- 

415 Jones fluid from the NpT plus test particle method. Mol. Phys. 76 (1992) 

416 1319-1333. 

417 [29] J.E. Mayer. The statistical mechanics of condensing systems. I. J. Chem. 

418 Phys., 5 (1937) 67-73. 

419 [30] J.E. Mayer. Statistical mechanics of condensing systems V Two-component 

420 systems. J. Phys. Chem., 43 (1939) 71-95. 

421 [31] J.E. Mayer and M.G. Mayer. Statistical Mechanics, John Wiley and Sons, 

422 New York, 1940. 

423 [32] C. G. Gray, K. E. Gubbins. Theory of molecular fluids, 1. Fundamentals, 

424 Clarendon Press, Oxford, 1984. 

425 [33] R.D. Mountain. A polarizable model for ethylene oxide. J. Phys. Chem. 

426 B, 109 (2005) 13352-13355. 

427 [34] M. p. Allen, D. J. Tildesley. Computer simulations of liquids. Clarendon 

428 Press, Oxford, 1987. 

429 [35] M. Ricci, M. Nardone, F. Ricci, C. Andreani, and A. Soper. Microscopic 

430 structure of low temperature liquid ammonia: A neutron diff'raction exper- 

431 iment. J. Chem. Phys. 102 (1995) 7650-7655. 

432 [36] J. Vrabec and H. Hasse. Grand Equilibrium: vapour-liquid equilibria by a 

433 new molecular simulation method. Mol. Phys., 100 (2002) 3375-3383. 

434 [37] H. C. Andersen. Molecular dynamics simulations at constant pressure 

435 and/or temperature. J. Chem. Phys., 72 (1980) 2384-2393. 

436 [38] B. Widom. Some topics in the theory of fluids. J. Chem. Phys., 39 (1963) 

437 2808-2812. 



17 



438 [39] I. Nezbeda, J. Kolafa. A New Version of the Insertion Particle Method 

439 for Determining the Chemical Potential by Monte Carlo Simulation. Mol. 

440 Sim., 5 (1991) 391-403. 

441 [40] J. Vrabec, M. Kettler, and H. Hasse. Chemical potential of quadrupolar 

442 two-centre Lennard- Jones fluids by gradual insertion. Chem. Phys. Lett., 

443 356 (2002) 431-436. 

444 [41] H. Flyvbjerg, H. G. Petersen. Error estimates on averages of correlated 

445 data. J. Chem. Phys., 91 (1989) 461-466. 



18 



Table 1; Parameters of the new ammonia model. The electronic charge is 
e = 1.6021 • 10-19 C. 



Interaction 


X 


y 


z 


a 




9 


Site 


A 


A 


A 


A 


K 


e 


N 


0.0 


0.0 


0.0757 


3.376 


182.9 


-0.9993 


H(l) 


6.9347 


6.0 


-0.3164 






6.3331 


H(2) 


-0.4673 


6.8095 


-0.3164 






6.3331 


H(3) 


-0.4673 


-0.8095 


-0.3164 






6.3331 
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Table 2: Vapor- liquid equilibria of ammonia: simulation results using the new 
model (sim) compared to data from a reference quality equation of state [25] 
(eos) for vapor pressure, saturated densities and enthalpy of vaporization. The 
number in parentheses indicates the statistical uncertainty in the last digit. 



T 


Psim 


Pcos 


Pshn 


Pcos 


p"- 

r Sim 


Pcos 






K 


MPa 


MPa 


mol/l 


mol/l 


mol/l 


mol/l 


k J / mol 


k J / mol 


240 


0.12(1) 


0.102 


40.26(1) 


40.032 


0.066(5) 


0.0527 


24.11(1) 


23.31 


280 


0.60(2) 


0.551 


36.98(2) 


36.939 


0.280(8) 


0.257 


21.56(1) 


21.07 


315 


1.65(4) 


1.637 


33.76(3) 


33.848 


0.74 (1) 


0.744 


18.96(2) 


18.57 


345 


3.37(4) 


3.457 


30.45(4) 


30.688 


1.55 (1) 


1.624 


16.19(3) 


15.79 


363 


5.22(5) 


5.101 


28.17(6) 


28.368 


2.56 (2) 


2.544 


13.93(5) 


13.65 


375 


6.37(6) 


6.485 


26.18(7) 


26.502 


3.17 (3) 


3.459 


12.48(6) 


11.89 


385 


7.88(5) 


7.845 


24.05(9) 


24.608 


4.27 (5) 


4.554 


10.49(9) 


10.08 


395 


9.54(7) 


9.422 


20.9 (1) 


22.090 


5.66 (9) 


6.272 


8.1 (1) 


7.66 
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List of Figures 



1 Saturated densities of ammonia. Simulation results: • new model, 

■ model from Kristof et al., O first model; — reference EOS [23] . 
critical data from simulation: O new model, □ model from Kristof 

et al.; + experimental critical point [27] |22 

2 Vapor pressure of ammonia. Simulation results: • new model, 

■ model from Kristof et al., O first model; — reference EOS [23]. [23 

3 Enthalpy of vaporization of ammonia. Simulation results: • new 
model, ■ model from Kristof et al., O first model; — reference 
EOS [53] Q 

4 Relative deviations of vapor-liquid equilibrium properties between 
simulation and reference EOS [33] {6z = {zsim— Zcos) / Zcos)- • new 
model, ■ model from Kristof et al., O first model. Top: vapor 
pressure, center: saturated liquid density, bottom: enthalpy of 
vaporization [25 

5 Relative deviations for the density between simulation and ref- 
erence EOS 113] i^P = (Psim — Pcos)/Pcos) in the homogeneous 
region: O simulation data of new model, — vapor pressure curve. 
The size of the bubbles denotes the relative deviation as indicated 

in the plot [26 

6 Relative deviations for the enthalpy between simulation and ref- 
erence EOS [13] {Sh = {hsim — /icos)//ioos) in the homogeneous 
region: O simulation data of new model, — vapor pressure curve. 
The size of the bubbles denotes the relative deviation as indicated 

in the plot [27 

7 Second virial coefficient: • simulation data of new model, — 
reference EOS [23] m 

8 Partition functions of ammonia: — simulation data of new model, 

o experimental data ^35] |2S 
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Figure 2: Eckl et al. 
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Figure 4: Eckl et al. 



25 






1 1 1 

e 


© 


1 1 









O ; 











o 


O 


o 


o- 


O 


o 


o 




o 


© 


® 


® - 



600 



700 
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Figure 8: Eckl et al. 
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